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Abstract 

When particulate suspensions are sheared, perturbations in the shear flows around the rigid 
particles increase the local energy dissipation, so that the viscosity of the suspension is effectively 
higher than that of the solvent. For bulk (three-dimensional) fluids, understanding this viscos- 
ity enhancement is a classic problem in hydrodynamics that originated over a century ago with 
Einstein's study of a dilute suspension of spherical particles.^ In this paper, we investigate the anal- 
ogous problem of the effective viscosity of a suspension of disks embedded in a two-dimensional 
membrane or interface. Unlike the hydrodynamics of bulk fluids, low- Reynolds number membrane 
hydrodynamics is characterized by an inherent length scale generated by the coupling of the mem- 
brane to the bulk fluids that surround it. As a result, we find that the size of the particles in 
the suspension relative to this hydrodynamic length scale has a dramatic effect on the effective 
viscosity of the suspension. Our study also helps to elucidate the mathematical tools needed to 
solve the mixed boundary value problems that generically arise when considering the motion of 
rigid inclusions in fluid membranes. 
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I. INTRODUCTION 



The dynamics of particulate suspensions in a viscous fluid are central to a variety of funda- 
mental scientific questions in hydrodynamics, soft condensed matter, and biological physics. 
A rather common and useful simplification of these studies replaces this heterogeneous sys- 
tem with a coarse-grained homogeneous one that has modified physical parameters such as 
viscosity. In bulk (three-dimensional) suspensions, understanding this change in viscosity at 
a quantitative level has captured the interest of researchers in various disciplines for over a 
hundred years. Beginning with Einstein,-i^i^ the basic physical interpretation of this result 
emerged: Under an externally imposed shear, the fiuid in the absence of the suspension 
adopts a spatially uniform shear stress and dissipates energy per unit volume proportional 
to that stress. The coefficient of proportionality is the bulk fiuid viscosity rj^D. With the 
addition of the particulate suspension this uniform shear stress becomes incompatible with 
the fiow boundary conditions at the surfaces of the particles, leading to more complex fiows 
surrounding the particles. These additional fiows cause additional energy dissipation in the 
fiuid. Thus, the coarse-grained homogeneous fiuid must have a higher viscosity than the 
original fiuid. This effective viscosity 77^^ must depend on the particulate volume fraction: 
The larger the volume fraction, the more energy dissipated by the suspension and thus the 
higher the effective viscosity. For a dilute suspension of spherical particles of radius a and 
number density n, Einsteini found that, to leading order, the effective viscosity depends on 
the volume fraction 03d = ^na^n as 

(1) 

This result has been extended to non-rigid droplets in a fluid^ and even non-spherical geome- 
tries,~i^ where changes in the numerical prefactor are found. In all these cases, the results 
apply only at low particulate volume fractions. Experiments flnd that Eq. (jlj) holds for 
03D ^ O.Ol.-i^ Above these concentrations the hydrodynamic interactions between particles, 
which are neglected in these calculations, become important. At such volume fractions one 
must consider these effects, as well as the possibility that the imposed shear flow changes 
the microstructure of the suspensions.-'^'^ 

In contrast to this tremendous effort in exploring the effect of flnite particulate concentra- 
tions on the viscosity of three-dimensional suspensions, comparatively little is known about 
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the analogous problem for fluid membranes and interfaces. The problem of membrane hy- 
drodynamics is complicated by the interactions of the essentially two-dimensional viscous 
membrane with the surrounding three-dimensional solvents. Because of this coupling, in- 
plane fluid momentum in the membrane is not conserved: Around a moving point-like 
particle in the membrane, momentum transfers to the surrounding fluids over a length scale 
set by the ratio of the 2D membrane viscosity ?7m to the 3D solvent viscosity r^sD. This 
Saffman-Delhriick length io ~ ?7m/'73Er^'^ makes membrane hydrodynamics qualitatively 
distinct from the usual three-dimensional hydrodynamics of bulk liquids, since the latter 
theory has no analogous inherent length scale (in the limit of vanishing Reynolds number). 

The existence of an inherent length scale in membrane hydrodynamics has profound im- 
plications on the transport properties of membranes. The mobility /i of a particle in an 
overdamped bulk 3D fluid is always inversely proportional to its size a, yU ~ as long 

as the Reynolds number remains small. In contrast, the mobility of a membrane-bound 
inclusion exhibits two drastically different behaviors as the particle size is varied, depending 
on the ratio a/io- When the particle size is small compared to io, the flows in the membrane 
dissipate much more energy than those in the surrounding bulk fluids, and the mobility 
only has a weak logarithmic dependence on the particle size.— '^^'^'^ Conversely, when the 
particle size is large compared to io, the flows in the bulk dissipate more energy. Not sur- 
prisingly, this leads to a mobility that, like its three-dimensional counterpart, is inversely 
proportional to a,— although the numerical prefactor is different. This complex dependence 
of mobility on particle size, as well as the related complex distance-dependence of hydro- 
dynamic interactions,— '^'^'i^ have been strongly supported by several experiments.— i^^*^ 
However, we point out that recent work has suggested that protein transport in lipid bilayers 
is more subtle than suggested by the original Saffman-Delbriick model.— i^^i^I In spite of this 
subtlety, understanding the effective membrane viscosity remains important for studies of 
the diffusive properties of transmembrane proteins. 

In this article we examine the effect of a flnite but small concentration of membrane- 
bound particles on the effective membrane viscosity. In essence we wish to flnd a relation 
analogous to Eq. ([1]) expressing the dependence of the effective membrane viscosity on the 
area fraction of membrane inclusions. We have several different motivations to study this 
problem. First, there is the fundamental question of how the Saffman-Delbriick length 
enters the coefficient of the area fraction term in the membrane version of Eq. ([T]). As with 
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the mobility of a membrane-bound inclusion, that coefficient should be a function of the 
dimensionless ratio e = a/^o for particles of radius a. Secondly, the effect of particulate 
suspensions on membrane viscosity addresses important biophysical questions regarding the 
dynamics of proteins embedded in the plasma membrane of cells. It is now well known that 
cell membranes are crowded environments in which the diffusive transport of transmembrane 
proteins controls, for example, cell adhesion and cell-cell signaling. These problems in protein 
dynamics are currently the focus of much study.^^ While much attention has been paid to 
how immobile obstacles in the membrane can hinder diffusion,— i2i comparatively little has 
been paid to how a suspension of mobile particles can have a similar effect by increasing the 
viscosity of the membrane. 

Finally, this study provides a simple physical system in which to explore a class of complex 
mathematical problems known as dual integral equations. Such problems arise generically 
in systems involving the transport of rigid inclusions in fluid membranes. The fundamental 
mathematical difficulty presented in these systems is that their behavior is governed by the 
solution to a mixed boudary-value problem. Physically, this arises from the two distinct re- 
gions in the system: The fluid regions of the membrane and the solid regions of the particle 
interior. In the former region, a stress continuity condition applies; that is, the internal 
stresses caused by the flows in the membrane must be balanced by the external stresses 
exerted on the membrane by the surrounding solvents. In the interior of the inclusion, on 
the other hand, the particle's rigidity supplies arbitrary constraint stresses to ensure that 
the entire inclusion executes only rigid body motions. As a result, the boundary condition 
in this region becomes a "stick" velocity balance condition. Each of these boundary condi- 
tions is expressed as an integral equation, so that the system requires two separate integral 
equations to be satisfied simultaneously in two non-overlapping domains of the membrane. 
Dual integral equations such as these have certainly been studied before,— i^^i^i^i^ even 
in the context of membrane hydrodynamics,—"^ but we believe that the methods are not 
widely known. To that end, we include an explanatory Appendix recapitulating the basic 
mathematical tools needed to solve these dual integral equations as they arise in membrane 
hydrodynamics. By mastering these tools here, we open the possibility of exploring numerous 
related problems, including the lubrication forces between two approaching transmembrane 
proteins or lipid rafts. 

To briefly summarize our results, we find that, like the mobility of a membrane-bound 
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inclusion, the effective membrane viscosity divides naturally into two regimes, characterized 
by the value of the dimensionless parameter e = a/lo. When e <^ 1, we find that the effective 
viscosity does indeed behave in a manner analogous to the three-dimensional result Eq. ([T]); 
that is, we find r/J^f = r7m[l + /(ct/^o)0]) where the coefficient of the area density term is a 
function of e. When e ^ 1, however, most of the dissipation occurs in the bulk solvents, 
so it is more appropriate to consider the inclusions as shifting the viscosity of the three- 
dimensional viscosity of these solvents. Here, we find a result nearly identical to the original 
Einstein result Eq. ([1]). In addition, our calculation provides a solution for the effective 
viscosity for arbitrary e that interpolates between these regimes. The full solution relies 
on numerical integration and matrix inversion, but we provide a closed form approximate 
solution that is exact in both asymptotic limits and has at most a small (< 10%) error for 
intermediate values of a/i^. 

The remainder of the paper is organized as follows: In Section [Til we determine the 
shear flows and dynamic pressures around an isolated membrane inclusion. These flows 
are then used in Section IIIII to determine the effective viscosity of a dilute suspension of 
such inclusions. We conclude in Section IIVI with a summary of our results and a discussion 
of problems in the hydrodynamics of rigid membrane-bound inclusions which we plan to 
address in the future using the knowledge we have gained here. 

II. ISOLATED INCLUSION 

Consider a flat, two-dimensional membrane (located at z = 0) consisting of a distinct 
fluid of viscosity ?7m. The membrane is surrounded by bulk fluids above {z > 0) and below 
{z < 0) whose shear viscosities are //+ and respectively; see Figure [U^a). We assume that 
all three fluids are incompressible and that all flows occur at low Reynolds number. Thus, 
the membrane velocity field v™ must obey the 2D incompressible Stokes equation: 

r/mVK + a+ I^^Q - a-,|^^Q - d^p^ = 0, (2) 

V ■ v"^ = 0, (3) 
while the velocity fields of the bulk fluids obey the incompressible 3D Stokes equation: 

VV^ = -VP^, V-v^ = 0. (4) 
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FIG. 1: (a) Schematic illustration of a membrane (viscosity r]^) containing a suspension of disks 
(radius a) and surrounded by bulk fluids above {z > 0, viscosity r/+) and below {z < 0, viscosity 
r/_). (b) Detail of an isolated inclusion in the membrane, viewed from above. The origin of the 
in-plane Cartesian (x, y) and cylindrical (r, 9) coordinates is the center of the inclusion. 

Here, Pm and P± are the membrane and bulk fluid pressures respectively, and af^ = 
V± [^i'^j' + ^pt] is the bulk fluid stress tensor. 

In general, any membrane flow field can be decomposed into three linearly independent 
normal modes, which correspond to the compression, bending, and shearing of the mem- 
brane. The out-of-plane bending deformations are decoupled at linear order from the in- 
plane flows. Since the focus of this paper is the dissipation caused by the in-plane flows, we 
ignore all bending deformations. Furthermore, we eliminate the compression modes by our 
assumption of the incompressibility of the membrane. This assumption is generally appro- 
priate for lipid bilayers. Thus, the hydrodynamic flows in the membrane can be decomposed 
purely into shear modes; that is, any membrane fluid velocity fleld can be written as a linear 
superposition of these modes. It is known that pure shear flows in a flat membrane generate 
no pressure in the surrounding bulk fluids,'-^ so we set = 0. 

In order to calculate the effective viscosity of the membrane, we follow loosely the deriva- 
tion of Einstein's result for the effective viscosity of a dilute three-dimensional colloidal 



suspension given in Ref. Il3|. We impose a simple shear flow in the absence of the particulate 
suspension and calculate the resultant dissipative stress in the system. These flows in the 
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membrane and surrounding solvents act as a probe of the viscous dissipative processes in the 
system. We then add a single rigid particle to the membrane and calculate the consequent 
perturbation to the flow fields. Using these flows, we calculate the average stress tensor in 
a dilute suspension of such particles in the membrane. By examining the terms that arise 
from the in-plane dissipative flows, we extract the effective viscosity. 

The simplest membrane shear flows generate constant (i.e. position- independent) stresses. 
Thus, we choose the "unperturbed" membrane velocity v™ to be of the form 

Vo,i{x,y) = S^,Caf3Xf^, (5) 

where Cap is a traceless symmetric tensor. Throughout this paper we use Greek indices for the 
in-plane (2D) vector component x, y and Latin indices for the 3D vector components x, y, z; 
the delta function 6^^ projects the Latin indices onto the Greek indices. The symmetry of the 
tensor [c^/? = C/3q] excludes flows corresponding the rigid rotation of the entire membrane. 
Such flows generate no dissipative stresses in the membrane and therefore are unnecessary. 
The vanishing trace [caa = 0] enforces the incompressibility constraint Eq. ([3]). 

Given the velocity field Eq. ([5]), we need to determine the bulk fluid flows and mem- 
brane pressure p™. The bulk flows are governed by the incompressible Stokes equation, 
Eq. (jll). The boundary conditions are given by the usual "stick" boundary conditions at 
the surface of the membrane, Vq (x, y, 0) = v™(a;, y), as well as the the 2D Stokes equation, 
Eq. ([2]). It is straightforward to show that the shear flow in the membrane Eq. induces 
the same shear flows in the bulk fluids: 

VoA^^y^^) = ^LCapXf^. (6) 

Here and throughout the paper, we use the vector field \{x,y,z) to represent the velocity 
field throughout all space: 



v{x,y,z) = < 



V (x, y,z) z < 

v'^{x,y) z = 0- (7) 
v~^{x, y,z) z > 



Finally, the membrane pressure vanishes for the unperturbed flows, pj}^ = 0. This solution 
satisfies Eqs. ([2])-(jl]), as well as the stick boundary conditions, and is thus the unique solution 
for the velocity field at all points in the system. 



We now introduce an isolated membrane inclusion, a rigid disk of radius a and of negligible 
thickness, at the origin of our coordinate system; see Fig. [I](b). Its presence perturbs the 
flows in the system and introduces new boundary conditions not satisfied by the unperturbed 
flows given above. Due to the linearity of the Stokes equation, we can write the total fluid 
velocity as v = Vq + Vi; that is, Vi is the "perturbative" flow field caused by the inclusion. 
It is clear from the rotational symmetry of Eq. ([6]) that the disk remains at rest: 

v-(r, 6) = -Vo-(r, 6) = -rc^pn^ r < a, (8) 

where r, 6 are the radial and angular variables, respectively, in cylindrical coordinates [see 
Fig. [H^b)], and Ua = Xa/r is the in-plane unit vector. Furthermore, the perturbative flows 
must vanish far away from the disk: 

lim Vi(r, ^, 2;) = lim Vi(r, ^, z)=0. (9) 

r— >oo 2— utoo 

The final boundary condition is given by the 2D Stokes equation, Eq. (j2]), which holds 
everywhere outside of the disk r > a. We note that if we had included rotational flows in 
the unperturbed membrane flows (i.e. if we allowed Cq,/3 to have antisymmetric parts) then 
the disk would simply rotate rigidly with the fluid, thus generating no additional sources of 
dissipation. 

Since we have chosen the unperturbed membrane flows v™ to be entirely composed of 
shear modes, and these modes are linearly independent of the bending and compression 
modes of the membrane, the perturbative velocity field v™ must also consist solely of shear 
modes. As a result, it can be written as an antisymmetric derivative of a scalar field: 

vTJr,e) = e^pdpMr,0), (10) 

where €0/3 is the antisymmetric unit tensor. Indeed, v™ can also be written in this form, 
with the scalar field 

2 

Mr,0) = -j[c.,yCos29-c,,sm29] = -r^Q{9). (11) 

We know from the linearity of the Stokes equation and the azimuthal symmetry of the 
disk that the angular dependence of v™ is set by the boundary condition Eq. (jH]). Thus, the 
angular dependence of ■?/'i(r, 6*) must be identical to that of ipo{r,9): ipi{r,9) = p{r)Q{9). 
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Using separation of variables, it is straightforward to show that the incompressibility con- 
straint in Eq. (jlj) and the stick boundary condition at the membrane surface z = imply 
that the bulk fluid velocities f 6*, z) have the form 

t;J,(r, e, z) = 5ih^{z)e^^dp^,{r, 9). (12) 

Then the Stokes equation for the bulk velocities, Eq. (jl]), becomes 

p{r) " A±(2) ^"^"^^ 

From the boundary condition Eq. ([2]), it is clear that we should choose exponential decays 
for h{z). Then the solution to Eq. ( fT3l) is given by 



Mr,0)h^{z) = a'e{e) / -^B{q)Mqu)e-''^, (14) 

Jo ^ 

where u = r/a, ( = \z\ /a, and J2{qu) is the second order Bessel function of the first kind. 
The Bessel function of the second kind Y2{qu) is also a solution to Eq. ( |T3|) .but it fails to 
satisfy the requirement of finite fluid velocities at r = 0. 

The function B{q) is a modified Hankel transform of the function p{u), i.e. the radial 
dependence of the scalar field ipi. In general, the kernel of these transforms is the product of 
a Bessel function J^iqu) with q^ for arbitrary real numbers m and p. Using the orthogonality 
and completeness of the Bessel functions, it can be shown that there is a one-to-one mapping 
of the function p{u) defined on the half-line < u < oo and its modified Hankel transform 
B{q).^'^ At the moment, B{q) is an undetermined function. Its form is determined by the 
boundary conditions in the membrane, which are given below. 

It is straightforward to show using Eq. that 



v,4r,e,z) = -6ij^ dqB{q)e-'^i 



J I X 2J2(gM) , 
qu ' 



(15) 



In order to determine the function B{q), we need to enforce the boundary conditions 
in the membrane. These boundary conditions are integral equations for B{q), because the 
velocity Eq. (fTS!) is itself an integral equation. Inside the disk - that is, for -u < 1 - we impose 
the stick boundary condition Eq. ([8]). Since this condition must be satisfied for arbitrary 6, 
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we can see from Eq. ( IT5|) that we obtain two separate integral equations: 

POO 

/ dqq-^B{q)J2{qu) = u\ u<l, (16) 
Jo 

POO 

/ dqB{q)Js{qu) = 0, u < I. (17) 
Jo 

Outside of the disk, we have the stress balance condition Eq. By taking the anti- 

symmetric derivative eap.d^ of this equation, we can eliminate the membrane pressure. Then, 
using Eq. (IT^ . we obtain the final integral equation: 

f dqq^B{q)J2{qu)\l + -]=Q, M > 1, (18) 
Jo L eJ 

where e = f^^, with the Saffmann-Delbriick length = ^^^^ . The parameter e is the key 

control parameter for the hydrodynamics of membrane-bound inclusions. When e ^ 1, the 

flows in the membrane dissipate much more energy than the induced flows in the bulk fluids; 

conversely, when e -C 1, the dissipation occurs primarily in the bulk. 

Finally, we obtain the membrane pressure using Eqs. ([2]), f|T5l) . and f|T8|) : 

p"(r, 6) = ^Caf^rianp / dq B{q) {q + e) 
4 Jo 

X [quJiiqu) - 2J2{qu)] . (19) 

The integral Eqs. (fT6|) - (fT8|) completely determine the modified Hankel transform B{q), 
which in turn determines the fluid velocities and pressures everywhere in the membrane 
and bulk fluids. However, finding the solution to these integral equations is difficult. The 
difficulty arises from the fact that this is a mixed boundary value problem: Inside the disk 
(0 < M < 1), we have a Dirichlet boundary condition that sets the total membrane velocity to 
zero; Outside the disk, we have a Neumann boundary condition that imposes stress balance 
across the fluid membrane. As a result, the boundary conditions Eqs. ( lT6|) - f|T8l) ultimately 
reduce to a pair of dual integral equations. Specifically, we must find the transform B{q) 
that satisfies Eq. (|T6|) inside the disk and Eq. f|T8|) outside the disk; we show in Appendix |X] 
that Eq. (IT7|) is redundant, because it is automatically satisfied by the solution to the dual 
integral equations. 

By contrast, consider a problem in which the boundary condition is given by a single 
integral equation that is valid over the entire region < u < oo. In this case, the in- 
tegral boundary condition is easily inverted using the mutual orthogonality of the Bessel 
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functions.— This is analogous to the well-known inversion of the Fourier expansions of a 
function. Indeed, if the size of the inclusion is very small, a ^ 0, we can approximate it 
by a point-like particle and ignore the velocity balance conditions inside the disk, Eqs. (ITB]) 
and f[T7|) . This limit, which is used often in membrane hydrodynamics,— greatly 
simplifies the solution. In the present problem, though, the finite size of the inclusion is 
essential in determining the effective viscosity of a suspension since it controls the suspen- 
sion's area density. Furthermore, one of the major motivations of this study is to gain a 
better understanding of the mathematical difficulties encountered in solving dual integral 
equations. 

The mathematical tools necessary to solve these dual integral equations are derived in 
Ref. 0; we summarize the necessary results in Appendix \M Briefly, we need to flnd a way 
to reduce the two modifled Hankel transforms, Eqs. f|T6|) and f|T8|) . into a single modified 
Hankel transform valid over the entire half-line < m < oo. Once we have accomplished this 
task, we can invert the remaining transform using the inverse modified Hankel transform. 
In order to combine the dual integral equations, we need to transform Eqs. (fT6l) and (|T8ll 
using operators that act on these integral equations entirely within their respective regions of 
validity, < m < 1 and 1 < n < cxd. In addition, these integral operators must possess simple 
convolution properties with the modified Hankel transforms. Such operators are known as 
the Erdelyi-Kober operators.— 

In Appendix \^ we define the modified Hankel transforms and the Erdelyi-Kober oper- 
ators; we also present the relevant inversion and convolution properties of these operators. 
Using these properties, the dual integral Eqs. ( fT6l) and ( fT8l) can be reduced to a single in- 
tegral equation, Eq. (lAlSp . It is convenient to re- write this equation in terms of spherical 
Bessel functions jn{u) = a/tt/ {2u)Jn+i/2{u), so that 







(20) 




Using the addition theorem for spherical Bessel functions. 



oo 




(21) 



m=0 
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we find that the even m terms cancel, leaving 

oo 

[u + e) B{u) = 5^6„(e)j2n-i(M), (22) 

n=l 

where the coefficients are defined as 

16e 2 f°° 

6„(e) = 6n,i + - (4r2 - 1) / dz zB{z)j2n-iiz). (23) 

To solve Eq. (!22l) . we convert it into a matrix equation for the coefficients 6„ by multiplying 
it by uj2i-i{u)/{u + e) and integrating. Using the orthogonality of the spherical Bessel 
functions, 

POO 

dqj2n-i{q)j2i-i{q) = ^(Jf^i^^ (24) 

we find 

oo „ 

5^6„(e)A^„,; = -5,,i, (25) 
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where 



oo 



M^,^ I rf,^_H!!zi^)^^. (26) 







Eq. (125!) is a matrix equation for the coefficients 6„; given the matrix inverse Mi^, its 
solution is trivial: 

bnie) = ^^r,n- (27) 

Given the coefficients fe„ the function B{u) can be found from Eq. ( I22ll . Then the mem- 
brane velocity Eq. (fT5l) becomes 



N ( 



t;i,„(r,^,0) = -}_6„(e)<( - -c^pnpV{u-n,l) (28) 



+ 



1 



V(m; n, 3) 



V(M;n,m)= / — —j2n^i{q)Jm{qu). (29) 



where 



Finally we turn to the pressure Eq. (1191) . It is straightforward to show using Eq. (1221) 
that all of the integrals in the resultant expression vanish, except for the n = 1 term in the 
sum [see Eq. (1B6I) ]: 

Pi {r, 0) = &i(e). (30) 
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FIG. 2: Dimensionless velocities [see Eq. (j3ip ] (a) vqi (t>) vi, and (c) v around an isolated in- 
clusion of radius a, as a function of the dimensionless distances f , f • The colorfield (b, c) is the 
dimensionless membrane pressure p^. All fields are calculated using e = 1 and Cxx = Cxy = 0.1. 



We can see from Eq. ([5]) that the tensor has units of [sec] ^. Therefore, the appropriate 
dimensionless quantities are 



r _ r 

Ca/3 = V = -V, p"^ = pf, 

a Vm 



(31) 



where r is the characteristic time scale in Caf3- 

Figure [2] shows the dimensionless membrane velocity fields vq) vi, and v and a color 
plot of the dimensionless pressure field p™. In order to compute these fields, the integrals 
/iAn,i and V{u]n,m) must be computed numerically, and the matrix Ai must be inverted 
numerically. The details of this procedure are provided in Appendix [Bl Figure [2](a) shows 
the unperturbed velocity vq) which clearly does not respect the boundary conditions at the 
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surface of the inclusion. The perturbative velocity field vi shown in Fig. [2](b) accounts for 
these boundary conditions. We can see that the perturbative velocity inside the particle is 
equal and opposite to the unperturbed velocity, causing the total velocity to vanish there and 
thus respect the boundary condition Eq. ([H]), as shown in Fig.[2I^c). In addition, the insertion 
of the particle into the membrane gives rise to regions of positive membrane pressure where 
the perturbative velocity flows out of the inclusion; conversely, regions of negative membrane 
pressure arise where the perturbative velocity flows into the inclusion. For Figure [2], we have 
chosen an intermediate value of the dimensionless parameter e = 1, i.e. we have set a = Iq. 
For different values of e, the velocity and pressure fields look qualitatively similar to those 
in Fig. [21 since the boundary conditions at the surface of the inclusion must still be obeyed. 
However, if we increase the viscosity of the membrane while keeping the particle size constant 
- that is, if we decrease e - the gradients in the perturbative membrane velocity field vi 
are decreased, causing this velocity to persist farther away from the inclusion (not shown). 
In addition, the magnitude of the pressure field increases. Conversely, higher values of e 
lead to more localized perturbative velocity fields and smaller membrane pressures. We can 
understand this behavior in the following way: As mentioned above, viscous dissipation in 
the membrane dominates in the limit of small e. As a result, large gradients in the membrane 
velocity field are unsustainable, causing the perturbative velocity field at the surface of the 
inclusion, which is required by the boundary conditions, to persist farther away from that 
inclusion as e is decreased. 

III. EFFECTIVE MEMBRANE VISCOSITY 

Armed with the results of the previous section, we now turn to computing the effective 
viscosity of a dilute suspension of membrane-bound inclusions. As discussed above, we use 
the stress tensor to probe the dissipative processes in the system. The effective membrane 
description of the suspension implies a coarse-graining of the system over length scales much 
larger than the size of the inclusions. Thus, we compute the stress tensor averaged over the 
entire volume of the system Vtot, which can be written as 




(32) 
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Due to the cylindrical symmetry of the problem, we choose Vtot to be a cylinder whose height 
iftot and radius i?tot are large. This volume includes the interiors of the solid inclusions. 
Within these regions, the stress tensor 0"*°* is not simply the fluid stress tensor; rather, it is 
the solid stresses in the inclusion caused by the fluids flows that surround it. 

Up until this point, we have been treating the membrane as a strictly two-dimensional, 
flat surface. In this model, the stress within the membrane enters the integral above as a 
delta-function at the membrane surface z = 0. It is convenient to avoid such a singularity 
when calculating the average stress tensor. To do so, we use an equivalent three-dimensional 
model of the membrane for which the stress is continuous at all points. Specifically, a two- 
dimensional membrane with a viscosity rj^ and two-dimensional pressure is equivalent to 
a thin, three-dimensional fluid of thickness h, viscosity rj^/h, and bulk pressure p^/h, in the 
limit of a vanishing membrane thickness h 0.iiii2iJ^ A schematic illustration of this 3D 
model is shown in Fig. [31 Thus, we can compute the integrals in Eq. (!32|) using the three- 
dimensional model and then take the membrane limit /i — to recover the two-dimensional 
membrane considered in the previous Section. 

Consider first the membrane in the absence of the particulate suspension, with only the 
unperturbed flows vq present. From the results of Section [Til if is straightforward to see 
that the stress tensor takes the form 



1-j 



^ai^iMzKlS, viz) 



?7_ z < 

< z < h 
rjj^ z > h 



Vm 

h 



(33) 



Then the integral in Eq. f[32l) is given by 



(?7+ + ?7_) Cap + -z—r]raCal3 



tot 



(34) 



We now turn to the particulate suspension. In the 3D membrane model, each inclusion 
is a solid cylinder whose height h is equal to the membrane thickness; see Figure [31 We 
anticipate that the average stress tensor for the suspension will have the same form as 
Eq. (!34l) . with the membrane viscosity rj^a being replaced by an effective membrane viscosity 

. Thus, we define the effective membrane viscosity via the average stress tensor: 

2 



—tot 



eff 



tot 



(35) 
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In order to calculate the average stress tensor for the suspension, we need the total velocity 
v*°* for this system. Since we work in the dilute limit, we can ignore the hydrodynamic 
interactions between the particles in the suspension. That is, we discard the negligible 
alterations of the flow fields around one disk due to the presence of the other disks in the 
suspension, so that each disk is treated in isolation. In this limit, the total fluid velocity is 
simply a linear superposition of the unperturbed membrane flows vq and the perturbative 
flows from each disk in the suspension: 

N 

^,'"\r,e,z)=Mr,0,z) + Y,^^'\r,e,z), (36) 

n=l 

where v^\r, 6, z) is the perturbation to the flows vq caused by an isolated disk whose center 
is located in the membrane at position x*^"-*, which can be obtained from Eq. (fT5|) by a simple 
coordinate translation. 

Rather than attempting to directly compute the average stress Eq. ( 132|) for the suspension, 
we first separate out the contributions of the unperturbed flows and of the perturbative flows 
of each particle in the suspension. This can be accomplished by writing the average stress 
tensor as 



where 

1 r „ 



a\f = r]{z)d,v^ + r]{z)djVi + H^,, (37) 



^"-^ ~ V , 



a)f-v{z){dd-' + d,vT') . (38) 



Consider the first two terms in Eq. fl?r|) . Clearly, the contributions of the unperturbed 
flows to these terms will yield the unperturbed average stress tensor ao,ij, Eq. ( l34l) . Fur- 
thermore, we can show that the perturbative flows v^"'' do not contribute to these terms. 
Specifically, consider the quantity ri{z)diV^j . This clearly vanishes for j = z, but it also 
vanishes for i = z, because angular integration over an odd number of in-plane unit vectors 
n will vanish. The i = a terms also vanish, because the integral evaluates to the velocity at 
Xa = ±oo, where it vanishes. Thus, the flows v^"^ do not contribute to the first two terms 
in Eq. §7^: 

^If = ^o,ij + ^ij, (39) 

We now turn to the integral Sj^. In the fluid regions of the system (i.e. outside of the rigid 
inclusions) the integrand is equal to the fluid pressure. However, we know from the results 
of the previous section that this pressure vanishes everywhere outside of the membrane. 

16 
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FIG. 3: Schematic illustration of an isolated inclusion in a thin layer of fluid of thickness /i, enclosed 
by a large cylinder of height H and radius R. The viscosity of the layer is rj^/h; In the limit /i — > 0, 
this system is equivalent to a two-dimensional membrane of viscosity ?7ni- The cylindrical volume 
V is divided into volumes above {V~^), below (V^), and within (14n) the fluid layer, as shown. The 
cylindrical caps on these volumes are denoted by and Cj^, while the shells are denoted by 
and Sm, as shown. 



Furthermore, in the fluid regions of the membrane, we see from Eq. f l30l) that the angular 
dependence of the membrane pressure ~ CapnaUp. Averaging over the angular variable 9 
produces the integral 



TT 



2tt 



(40) 



From this we see that the contribution of the fluid membrane regions to Sj^ also vanishes 
since Cajj is traceless. Thus, the only regions of integration that contribute to Sjj are the 
solid interiors of the disks themselves. Due to our neglect of the hydrodynamic interactions 
between the disks (as justified by the assumption of a dilute suspension), each disk in the 
membrane provides an identical contribution to Hj^, so we have 

N 



-f^tot^tot JV 



(TX 



(41) 



where Atot = T^Rtot ^ is the number of particles in the suspension. The perturbative 
stress tensor ai^ij = aij — cro,ij, where, aij is the stress tensor everywhere within a system 
containing an isolated inclusion. Thus, we have converted the computation of the average 
stress tensor of a particulate suspension into the problem of a single isolated inclusion 
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considered in Section [TTl Although the integrand is non-zero only within that inclusion, it 
proves useful to re-extend the region of integration V to include all of the surrounding fluids. 
Therefore, we choose to be a large cylinder whose height H and radius R will eventually 
be taken to infinity; see Fig [31 

Consider the first term in Eq. fl4ip . the integral of the perturbative stress tensor ai ^j. 
From its definition, we can see that ai^ij contains all of the solid stresses within the inclusion, 
as well as the fiuid stresses caused by the perturbative velocity field vi. Since the integration 
domain V in Eq. (14T|) clearly includes the interior of the solid inclusion, we would need to 
determine the solid stresses in this region to compute this integral directly. We can avoid 
this difficulty, however, by using the divergence theorem to convert this volumetric integral 
into a surface integral. Stress continuity requires that dkCTik = at all points in space, 
including the interior of the inclusion. Furthermore, it is clear from Eq. fl33|) that dkcro,ik = 
everywhere. Then we may write 



Here, we have extended the height H of the enclosing cylinder to infinity. Because of the 
exponential decay of the perturbative fiuid velocity Eq. (IT^ as z ±oo, we neglect the 
integration over its circular end-caps at z = ±H/2 (see Fig. [3]). 

For i = a,j = z or i = z,j = a, it is straightforward to show - using Eq. f|T5l) and 
the fact that (9^^1,7 is even in n - that the integrand of the surface integral in Eq. (142|) is 
odd in h, and therefore vanishes upon integration over 6. For i = j = z, the integrand is 
oc c-ypri-ynp, which also vanishes upon integration over 9, by Eq. (HUj) . Thus 



The remaining terms in Sj^ are proportional to the discontinuous viscosity function ri{z). 
For these terms, we break up the integration volume V into three different regions containing 
the three separate fluids in the system. Namely, we divide V into three separate cylinders 
V~^,V~, and V™, which enclose the regions z > h, < z < h, and 2; < 0, respectively; 
see Fig. |3l Using the divergence theorem, we obtain integrals of the velocity components 




(42) 




(43) 
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Vi,a over the cylindrical shells S~^, S~ , and 5*™, whose outward normals are all n, as well as 
integrals over the end-caps Cj^, whose outward normals are cx z. Because the velocity is odd 
in h [see Eq. (IT51) ]. the latter integrals will all vanish. Thus, 



(44) 



V 



dzrj{z) 



2lT 



r=R 



From Eqs. (fT5|) and (fT9l) . we see that the z integrals in Eqs. (H3il and 
Returning to the limit of an arbitrarily thin membrane, /i — > 0, we find 

e 



lim 



ri{z)e 



-q\z\/a 



1 + 



are all identical. 



(45) 



From Eqs. (H3l) - P5|) . we find that Eij cx S^^S^j/Htot. Thus, the average stress tensor 
Eq. ( !39l) does indeed take the form of Eq. ( !35l) . as anticipated. Specifically, if we compute 
the remaining angular integrals in Eqs. (H3i) and (l44l) using Eqs. (ITSl) and (|T9l) . we find that 
the effective viscosity is 



eff 



Vn 



1 + 6 lim $ 



R 



where = A^vra^/Atot is the area fraction of particles in the membrane and 



^(u) = u 



dqB{q){q + e) 



(46) 



(47) 



In Eq. fl46|) . we have taken the radius R of the enclosing cylinder V to infinity, as promised. 

We have succeeded in finding an expression for the effective membrane viscosity in terms 
of the function B{q) determined in Section [Til Using Eq. fl22|) . we may write Eq. fj471) as 



n=l 



dqi2n-i{q) 



(48) 



It is straightforward to show [see Eq. flB6P ] that the integral of the second term in brackets 
vanishes for all ri, while the first integral vanishes for all n > 1. The n = 1 integral is cx I/m^, 
so its contribution to $(«) is -u-independent and thus survives the R ^ oo limit in Eq. ( I46l) . 
Thus, the effective membrane viscosity is determined solely from the coefficient &i(e): 



r]^ [1 + 0/(e)] . 



(49) 



The function /(e) can be computed numerically for arbitrary values of e; see Appendix [Bl 
Before we discuss this solution, however, we first consider the asymptotic limits e ^ 1 and 
e ^ 1, where analytic solutions for /(e) can be obtained. 
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A. Iq: Large Inclusions 



When the particle size a is much larger than the Saffman-Delbriick length io - that is, 
when e ^ 1 - viscous dissipation occurs predominately in the surrounding 3D fluids, rather 
than in the membrane. In this limit, we calculate the leading order and next-to- leading 
order dependence of 61(e) on e. We write 



The leading order term bn'' is found by approximating u + e ^ e in Eq. 
orthogonality of the spherical Bessel functions Eq. (I24l) . Eq. (I25l) becomes 



n=l 



2e(4/ - 1; 



16e 



TT 



ki- 



rn 

Using the 
(51) 



Using Eq. (!50i) . the next-to- leading order terms in Eq. (!25ll are 



(52) 



n=l 



where bn = bn/e, ^ = 1/e, Ai^l given by Eq. fl2^ . and 

poo 



A^^'J(0 = lim 



udu . ( ■ / 

^ I ^~ J2n-1 7 J2I-I 7 



(53) 



The region < m < ^ of this integral gives a negligible contribution in the limit ^ — > 
and can be discarded. In the remaining integral, we expand the spherical Bessel functions 
for large values of their arguments: j2n-i{x) ~ (— l)"cos(a;)/x for a; ^ 1. Using this 
approximation we find that the dominant contribution is logarithmic: 



1 



\n+l 



(54) 



Eq. (!52|) must be satisfied for arbitrary values of ^, as long as ^ is sufficiently small. As a 
result, the two terms in the sum of Eq. fl52|) must have the same functional dependence on 
C,. Then we find 

^«(0 = ^^^^(-l)"ein(0, (55) 

Thus, in this limit 

/(^) = ^ + ^Me), (56) 



TT TT^ 
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As mentioned above, the flows in the bulk fluids dissipate much more energy than the 
flows in the membrane in the limit e ^ 1. Therefore, it is more appropriate to deflne an 
effective three-dimensional viscosity in this limit. If we return to the average stress Eq. 
for a symmetric membrane (?7_|. = ''73d); 

2io 8a(j) " 



1 + 



(57) 



In the limit 0, the second term vanishes. If we compare this to the average stress of 



the unperturbed membrane Eq. (iMll in this limit, 

aij = 2S^iS^jri3D, (58) 
we can see that the effective 3-D viscosity is 



^730 = ^730 



1 + -03D 



(59) 



where 03d = Nna^ /Vtot acts as a volume fraction of the particles. This identiflcation is not 
precise, since we have taken the z-extent of both the membrane and its inclusions to vanish 
in order to compute the effective membrane viscosity. In writing (p^D in this form, however, 
we have given the inclusions a vertical (z) size of a and neglected any numerical prefactors 
of order unity associated with the precise geometry of the inclusions (e.g. cylinders vs. 
spheres). It is interesting to note that, in spite of this imprecision, the numerical prefactor 
S/tt ^ 2.55 is within 2% of the Einstein coefficient of 5/2 for a three-dimensional suspension 
of spheres, Eq. ([1]). This result is not unexpected, despite the fact that there are profound 
differences in the underlying assumptions regarding the particle distribution between our 
calculation and that of Einstein. In our model, all of the particles are confined to a plane, 
while in Einstein's work, the particles are assumed to occupy all space. This distinction is 
lost, however, in taking the low volume fraction limit, where both calculations reduce to a 
single-particle calculation. In this essentially mean- field limit, all hydrodynamic interactions 
between the particles are ignored, and the effective viscosity can depend on the mean volume 
fraction of the particles alone and not the details of their spatial distribution. Thus, it is 
not surprising that our result closely approximates Einstein's when the dissipation in the 
membrane is negligible. 
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B. a <^ £q: Small Inclusions 



We now consider the limit in which the Saffman-Delbriick length is large compared to 
the size of the inclusion, a <^ Iq. We still assume that the suspension is dilute, so that is 
small compared to the mean lateral separation of the inclusions, io <^ aj^f^- In this limit, 
we may continue to neglect the hydrodynamic interactions between the inclusions. We now 
expect that the viscous dissipation occurs predominantly in the membrane. 

In this limit, we expand 

^^--4, 6„(e)^6W+e6W. (60) 
We require that Eq. (12^ be satisfied term by term in e, so that 

= Ir^. (61) 



6(1) = 3 E RT^Ri'nKi^ (62) 

l,n=l 

where the matrix elements Rn^i and Rnj are given by Eqs. 
matrix inverse ^le computed analytically; see Ref. 

by 



B2l) and (lB3p . respectively. The 
l4 The elements i?|f ^ are given 



vrn \ r (n 

Then 61°) = 12 and 



Thus, in the limit a <ti io (i-e. e ^ 1), we find 



5e 



/(e) = 3 + -, e«l, (65) 



so that, using Eq. (I49l) . 



TT 



1 + 30 + —. 

7r£o 



(66) 



The leading-order correction to the viscosity, which is reminiscent of the Einstein result 
Eq. ([T]), gives the dependence of the effective membrane viscosity on the area fraction in 
the limit of vanishing inclusion radius. The next-to-leading order contribution shows an 
additional viscosity enhancement, proportional to a/io. This term is due to the additional 
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0.01 0.1 1 10 50 

FIG. 4: Exact numerical results (dots), asymptotic limits (dashed lines), and interpolation function 
(solid line) for the function /(e). For the numerical results, the infinite matrix is truncated at 6 x 6, 
and Q = 5e; see Appendix [Bl 

dissipation caused by the flows induced in the bulk fluids for small (but flnite) values of 
a/4. 

It is important to recognize the distinction between the e = limit of our system and 
the strictly two-dimensional system of a suspension of infinite cylinders in a bulk fluid that 
has been considered previously.— Even when the viscosities of the surrounding bulk fluids 
vanish (i.e. £o = oo), our system is not two-dimensional: The flows in the membrane extend 
inflnitely far away from the membrane in the bulk fluids, but the pressure is still non-zero 
only within the membrane. In contrast, the fluid pressure around an inflnite cylinder is the 
same everywhere along the cylinder axis. Thus, the e = limit of our system, where the 
viscosity correction = 30, is different than the viscosity correction of 20 for a suspension of 
cylinders.— 1^ 

C. Arbitrary a/lQ 

While the above asymptotic results are suggestive, it is clearly desirable to examine /(e) 
for arbitrary values of e. Figure H] shows the function /(e) over several decades of e values. 
The exact numerical solution to Eq. (H9l) is indicated by the points (see Appendix [B] for 
details), while the asymptotic values Eqs. fl56l) and fl65l) are indicated by the dashed lines. It 
is clear that the numerical solution agrees with the asymptotic expressions in the appropriate 
limits, and that the transition between these two limits is smooth and monotonic. 
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We can use the analytic expressions for the small and large e behavior of /(e) to construct 
an analytic function /(e) that smoothly interpolates between these extremes. We note from 
Eqs. dSSD and (jHSD that /(e) grows as e + ln(e) for large e but has no logarithmic divergence at 
small e. Thus, we need a term in the interpolation function /(e) that grows logarithmically 
at large e but remains well-behaved (i.e. non-singular) at small e. The obvious choice is 
ln(l + e): 



/(e) = ^ln(l + e)+^7(e), (67) 



where 



3+(^-^)e e«l 
g{e)={ . (68) 

e > 1 



To find a suitable function g{e), we employ the method of two-point Pade approximants:— 

= (69) 

where Ajq and Cm are polynomials of order N,M, respectively, in e. Without loss of 
generality, we can set Ca/(0) = 1. Therefore, we have N+M+1 unknown coefficients. Three 
of these coefficients can be set by the known asymptotic limits of g{e) given in Eq. fl^Sl) . 
Furthermore, in order to obtain g{e) ~ e for e ^ 1, we must have N = M + 1. Thus, the 
first non-trivial Pade approximant for g{e) is N = 2, M = 1: 

ao + aie + age^ 

This has four unknown coefficients, so there is not a unique Pade approximant for this 
function. However, it is straightforward to show that Ci is the undetermined coefficient, and 
that all values of Ci > give a smooth, monotonically increasing Pade approximant, so we 
set Ci = 1 for simplicity. By expanding g{e) for large and small values of e and matching 
these limits to those given in Eq. ( l68l) . we find 



/(e) = - In (1 + e) + ^ . 71 

Vr^ TT^ (i -|- €) 

In Figure m we plot the interpolation function /(e) as a solid line. We see that it exhibits 
excellent agreement with the exact numerical results for all values of e. Indeed, the error 
between /(e) and the numerical solution never exceeds 8%, as shown in the inset of Fig. HI 
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IV. CONCLUSION/FUTURE WORK 



Membranes and fluid interfaces are by tlieir nature liybrid systems. Althougli tlie mem- 
brane/interface is itself two-dimensional, it is surrounded by bulk three-dimensional fluids. 
As a result, the hydrodynamics of membranes can exhibit both two-dimensional and three- 
dimensional characteristics, depending on the system in question. The effective viscosity of 
a fluid membrane containing rigid inclusions demonstrates this dimensional crossover. For 
particles whose radii a are small compared to the Saffman-Delbriick length io, the effect of 
the suspension on the large length scale viscous dissipation under shear can best be thought 
of as providing an increase in the effective viscosity of the membrane rj'^. Conversely, for 
large inclusions relative to Eq, their effect on the sheared membrane and surrounding solvent 
can be understood as an increase in the viscosities of the bulk solvents that is proportional to 
the volume fraction of the inclusion. For arbitrary inclusion size, we have determined a rea- 
sonably simple interpolation formula that gives an accurate estimate of the exact numerical 
solution for the effective membrane viscosity. 

Mathematically, we have seen that the hydrodynamics of membranes containing rigid 
inclusions is a mixed boundary value problem whose solution obeys a set of dual integral 
equations. One of the benefits of this work is that it helps to elucidate the mathematically 
machinery needed to solve these dual integral equations. Since these equations arise in many 
problems in membrane hydrodynamics, we plan to use our newfound understanding of the 
mathematics to solve other problems. Despite the mathematical complexity of the problem 
presented here, it is in a way one of the simpler problems one can consider in the hydrody- 
namics of membranes with rigid inclusions, because the hydrodynamic interactions between 
the inclusions can be ignored. Indeed, one of the motivations for this work was to study a 
system in which we could learn about membrane hydrodynamics with rigid inclusions with- 
out the additional complication of particle interactions. Armed with this knowledge, we 
plan to investigate such interparticle interactions in future studies. One problem of partic- 
ular biophysical significance is the study of lubrication forces between two large membrane 
inclusions in close proximity. It is now widely believed that many transmembrane proteins 
recruit lipid rafts^ in the cell's plasma membrane. Treating these extended structures as 
essentially rigid objects, one may ask how the hydrodynamic interactions between two such 
rafts affect the kinetics of protein aggregation in the membrane. Similar questions can also 
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be asked of the kinetics of phase separation in the multicomponent hpid bilayers of giant 
unilamellar vesicles.— More generally we expect the mathematical and physical features of 
the problem considered in this paper to arise in the study of the kinetics of inclusions or 
finite size domains in any lipid bilayer or Langmuir monolayer system. Such problems should 
exhibit the phenomenon of a scale-dependent dimensional crossover explored here. 

MLH and AJL thank H. Stone for interesting conversations. This work was supported in 
part by grant Insf-cmmi 0800533t 



APPENDIX A: DUAL INTEGRAL EQUATIONS 

In this Appendix, we present the mathematical tools necessary to manipulate the integral 
equations, Eqs. f|T6|) - f|T8|) . and then use these tools to derive a single integral equation. The 
necessary integral operator identities are presented here without proof; we refer the reader 



to Ref. 
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for the derivation of these identities. 



Consider a function f{q) defined everywhere on the positive-g axis, < g < oo. Adopting 



29, we denote the modified Hankel transform of this 



the compact notation used in Ref. 
function by the operator Sr^^\f{u), which is defined by 

s,,./(«)^i;rd,%M;(,). (Al) 

"Jo y 

Where necessary we use the expanded notation Sr^^xf{u) = Sr^^x {fil)] u}. Using this nota- 
tion, Eqs. (fT6|) and (fTSjl can be written as, respectively, 

So,2Biu) =4, M < 1, (A2) 

Ss.B(u) = -—S2-2B(u), u>l. (A3) 

2' ex ' 

We will return to the final integral equation, Eq. (ITTI) . at the end of this Appendix. 

The principal difficulty presented by Eqs. ( ]A2[) and ( 1A3I) is that the unknown function 
B{u) is defined by two separate integral equations, each with its own domain of applicability. 
From the closure relation for Bessel functions, we know that that the inverse of a modified 
Hankel transform is another modified Hankel transform: specifically, S~\ = Sn+x-x- There- 
fore, the inversion of a modified Hankel transform is possible only if it appears in an equa- 
tion that applies to the entire half- line < u < oo. The dual integral equations, Eqs. flA2p 
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and (lA3p . clearly do not satisfy this requirement. As a result, we cannot directly invert the 
modified Hankel transforms in Eqs. flA2l) and flA3l) to solve for the function B{u). 

To resolve this dilemma, we combine Eqs. (lA2p and (1A3P into a single integral equation 
using the Erdelyi-Kober operators defined below. Through the application of these opera- 
tors, we can write the left-hand side of Eqs. (lA2p and (lASp in the same form. In this way, 
we generate a single integral equation whose domain of validity extends over the entire real 
positive axis. 

The Erdelyi-Kober operators are defined as 



Kr,,xf{q) 



2q 



2r] 



u 



r(A) 

These integrals only converge for A > —1/2; for A < —1/2, 

Ir,,\+nf{q) 



f{u) du. 



where n is an integer such that X + n > and 



-2r;-2A-2n-l 



q^ 



K, 



rj—n,\+n 



f{q) 



q 



2rj-2n-l 



^^^^^^ = 2d-q 



(A4) 
(A5) 

(A6) 
(A7) 

(A8) 



The utility of these operators stems from the following observations: (i) When acting 
on a function /(g), the operators /^,a and Kri^x involve integrals over (0,g) and {q,oo), 
respectively. Thus, they depend on two disjoint subspaces of the positive real line. This is 
essential because it allows one to apply In^x to Eq. ( ]A2p and obtain an integral that is well- 
defined for < M < 1. Similarly, we can apply \ to Eq. (]A3P and obtain an integral that 
is well-defined for 1 < m < oo. (ii) Both operators /^,a and Kn,\ have the simple convolution 
properties with modified Hankel transforms, namely 



(A9) 



and 



K„\S, 



ri,X'Jrj+X,'y — ^r;,A+7 



(AlO) 
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We now apply / to Eq. (lA2p and K to Eq. (lA3p . making a judicious choice of the coeffi- 
cients so that the left-hand sides of the resulting equations are identical. Specifically, it is 
straightforward to show using Eqs. (1A9P and (lAlOp that 

l2-^So^2 = KqzSs = SqI. (All) 

The coeffients of these operators are unique and set by the coefficients in Eqs. ( 1A2I) and (lASp . 
along with the requirement that the resultant modified Hankel transform be the same in 
both equations. Thus, Eqs. (1A2I) and (1A3I) can be written as a single integral equation. 



S,.B{u) 



-fi^O,! [^52,-2i?H] U>1 

From Eqs. (IA4p and flA6p . we can show that 



(A12) 



^2.-1 (4) 



^'^2 1 (4) 



150F 



16 



Using the identity 



along with Eq. ( lAlOp . we find 



^0^ 

u, 2 



-52,-25(^i) 



u 



-Ki ■sS2^2B(u) = -Si _iB(u) 

M 2-2 ' 2' 2 ^ ' 



(A13) 



(A14) 



(A15) 



Finally, by applying the inverse Hankel transform S \ = Si _i, the single integral equation 

U, 2 2 ' 2 

reduces to 

w <1 



Biu) = Si _i 

2 ' 2 
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U 



(A16) 



-5i i^C^x;) > 1 



We can write Eq. flA16P as a conventional integral equation using Eq. flAip . Noting that 
Ji/2{wu) = a/2/ ijiwu) sin(wM) , we find 



TT 



sm M cos M 



u 

2 



(A17) 



en 







dz zB{z) / (itu sin(wM) sin(wz). 
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Writing sin(wu) sin{wz) = cos[w{u — z)] — cos[w{u + z)] and noting that dw cos{wy) 
TT6{y), we perform the integration over w, 

sin u cos u 



71 (u + e)B(u) = 16e 



dz zB{z) 



u 



sinlM — z) sinfw + z) 



(A18) 



u — z 



u + z 



It is this form of the single integral equation that we use in Section [TTl to solve for the function 
B{u). 

Finally, we return to the third boundary condition, Eq. f[T7|) . which to this point we have 
neglected. For this integral equation the kernel is Js{qu). It is straightforward to show that 

j[°°drtw(5^-^)=0, (A19) 

Similarly, for all w > 1, 

poo 

dqJsiqu) sin(wg) = 0. (A20) 



From these two integrals we note that Eq. f[T7j) is automatically satisfied by any solution 
B{u) of Eq. (]A17|) : that is, this condition provides no unique information about the function 
B{u). 



APPENDIX B: NUMERICS 



In order to obtain a numerical solution for the coefficients 6ri(e), we need to compute 
the integrals A4n,i defined in Eq. fl2B]) . However, the oscillatory nature of the spherical 
Bessel functions makes these integrals difficult to compute numerically. We can avoid this 
difficult by dividing the integration region into two portions: < q < Q and Q < q < oo, 
where Q ^ e. Then in the latter region we perform a Taylor expansion of the denominator. 
Keeping terms up to second order in 1/q, 



Mn,l 



Q 



dqj2n-i{q)j2i~-i{q) 



+ Rn,l — ^Rn,i 



1 



q 



- + - 

q r. 



where^ 



Rn,i = [ —32n~i{.q)32i~i{.q) 

Jo q 



(Bl) 



(B2) 



2 (4 (n - ly -l){n + l-l){n + l) 



29 



f°° dq 

Rn,l = / — j2n-l(g)j2Z-l(g) (B3) 

Jo 1 

= -a{l + i) 6n,i+i + a{l)6n,i + |a (/ - |) Sn,i--i, 

with 

"^^^ " (4/ + l)(4/-l)(4/-3) - 
The remaining integral in A4n,i, whose integration region is < q < Q, can be computed 
numerically, as long as e is not too large (e < 50). This procedure yields accurate values for 
these integrals for sufficiently large values of Q; throughout this paper, we use Q = lOe. 
To compute the membrane velocity Eq. fl28|) . we proceed in the same manner: 
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1 1 e 
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g + e q q^ 



+ 2{u;n,m,l) — el{u;n,m,2), (B5) 

where for m > 1 and — | < p < m + 2n, 

/■°° dq 

I{u;n,m,p) = J — j2n-i{q)Jm{qu) (B6) 

_ V^r{n + ^)F{n+^,n-rn±P-2n + l-^) 
2P+iM2n-pr (2n + i) r - n + 1) 

and F{a, (5]^] z) is the hypergeometric function. It is very important to note that this 
integral vanishes if |(m + p) — n + 1 is a non-positive integer, due to the divergence of the 
Gamma function in the denominator. 

Finally, in order to determine the coefficients 6n(e), the infinite matrix Ain,i must be 
inverted. To do this inversion numerically, we truncate the matrix. It is straightforward 
to verify that accurate solutions for the velocities and effective viscosity are obtained for 
reasonably-sized matrices; for all of the numerical results presented in this paper, we truncate 
the matrix at 10 x 10. 
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